/*********** NOTES HEADER *********** 

DESCRIPTION: This program estimates the limited specifcation logit
	and estimates the marginal effects

Datasets used:
 (1) MatchedSample_Full

Datasets created:
 (1) Overall_diab_lim 
 (2) Overall_hyp_lim 
 (3) Overall_death_lim 
 
************************************/ 

cd "N:\MedicareClaims-P045601-BE\Work\hosp_retro\health_out\Program\PropensityScore\Patients\ResponsetoReferee\OtherPropScoreMethods\BiasAdjustment"
set more off
set matsize 10000
capture log close
clear all

global origData "N:\MedicareClaims-P045601-BE"
global dataIn "N:\MedicareClaims-P045601-BE\Work\hosp_retro\health_out\Data-In\"
global dataOut "N:\MedicareClaims-P045601-BE\Work\hosp_retro\health_out\Data-Out\"
global dpath "N:\MedicareClaims-P045601-BE\Work\ay_data"
global dataProp "N:\MedicareClaims-P045601-BE\Work\hosp_retro\health_out\Data-Out\PropScore\Patients"
global logs "N:\MedicareClaims-P045601-BE\Work\hosp_retro\health_out\Logs\PropScore\Patients"
global dpath "N:\MedicareClaims-P045601-BE\Work\ay_data"
global skapath "N:\MedicareClaims-P045601-BE\Work\ska"

adopath +  "N:/SIL-Common/estout"
adopath +  "N:/SIL-Common/outreg2"
adopath +  "N:/SIL-Common/reghdfe-master/package"

log using "./1.PatScore_overall_noweights.log", replace

	program getColnames 
	 quietly {
	  matrix A = e(b)'
	 
	 local rownames : rowfullnames A
	 local c : word count `rownames'
	 local names : colfullnames A

	 svmat A, names(col)

	 gen rownames = ""
	 forvalues i = 1/`c' {

	  replace rownames = "`:word `i' of `rownames''" in `i' 
	   } 

	 }

	end


	program getRegStats
	 quietly {
	  matrix A = e(b)'
	 
	 local rownames : rowfullnames A
	 local c : word count `rownames'
	 local n = `c' + 1
	 local pll = `c' + 2
	 local rsq = `c' + 3
	 local ndrop = `c' + 4

	 replace y1 =e(N) in `n' 
	 replace y1 =e(ll) in `pll' 
	 replace y1 =e(r2_p) in `rsq' 
	 replace y1 =e(N_drop) in `ndrop' 

	 }

	end

use  ./MatchedSample_Full 
 
 
 ** START STATISTICS **

*Diabetes
	foreach var of varlist oc_diabcomp1 oc_diabcomp2 oc_glaucoma  {
	*This limits the sample to diabetes patients
	 forval i=3(1)3  {

	if regexm("`var'","oc_([a-z0-9]+)$") { 
	  local pre = "`=regexs(1)'"

	 capture drop zip3`pre' zipct`pre'
	 gen zip3`pre'=zip3dig 
	 egen zipct`pre'=count(zip3dig) if kv_`pre'==1 , by(zip3dig)
	 replace zip3`pre'=0 if zipct`pre'<1250
	 replace zip3`pre'=0 if kv_`pre'==0 
	 replace zip3`pre'=STATE_encode+1000 if zip3`pre'==0 
	 
	 capture rename mxicd3_chronic mxicd3_endocrine 
	 
	 logit `var'  i.postmerger#i.sex_pat i.vmerger##i.sex_pat i.vmerger#c.nyear i.p_age i.agecat##i.sex_pat i.qtr [pweight=samp_weight] if flag`i'==1 & kv_`pre'==1 , cluster( propScoreId)
	 
	   matrix sesq = vecdiag(e(V))' 
	  svmat sesq
	  getColnames
	  getRegStats

	  rename rownames rn`pre'_`i'
	  capture drop j 
	  capture drop i
	  rename y1 beta_`pre'_`i'
	  rename sesq1 sesq_`pre'_`i'

	  nlcom exp(_cons+ _b[1.postmerger#1.sex_pat]+_b[1.vmerger] + _b[77.p_age])/(1+exp(_cons+ _b[1.postmerger#1.sex_pat]+_b[1.vmerger] + _b[77.p_age])) - (exp(_cons+  _b[77.p_age]+_b[1.vmerger])/(1+exp( _cons+ _b[1.vmerger] +_b[77.p_age])))
	   mat A=r(b)
	   mat B=r(V)
	   mat C=[A,B]
	  nlcom exp(_cons+ _b[1.postmerger#2.sex_pat]+_b[1.vmerger]+_b[1.vmerger#2.sex_pat] + _b[77.p_age] + _b[2.sex_pat]+ _b[3.agecat#2.sex_pat])/(1+exp(_cons+ _b[1.postmerger#2.sex_pat]+_b[1.vmerger]+_b[1.vmerger#2.sex_pat] + _b[2.sex_pat] + _b[77.p_age] + _b[3.agecat#2.sex_pat])) -   (exp( _cons+ _b[1.vmerger]+_b[1.vmerger#2.sex_pat]+_b[77.p_age] + _b[2.sex_pat] + _b[3.agecat#2.sex_pat])/(1+exp(_cons+ _b[1.vmerger]+_b[1.vmerger#2.sex_pat] + _b[2.sex_pat]+_b[77.p_age] + _b[3.agecat#2.sex_pat]))) 
	   mat A=r(b)
	   mat B=r(V)
	   mat D=[A,B]
	   mat E=[D',C']
	  svmat E
	  rename E1 mef_`pre'_`i'
	  rename E2 mem_`pre'_`i'

	   }
	  }
	 }

 save "./Overall_diab_lim.dta", replace

rename mxicd3_endocrine mxicd3_chronic 

*Hypertension
		foreach var of varlist oc_acuteicd8 oc_ami oc_ihd  {
		*This limits the sample to hyptertensive patients
		 forval i=2(1)2  {

		if regexm("`var'","oc_([a-z0-9]+)$") { 
		  local pre = "`=regexs(1)'"

		 capture drop zip3`pre' zipct`pre'
		 gen zip3`pre'=zip3dig 
		 egen zipct`pre'=count(zip3dig) if kv_`pre'==1 , by(zip3dig)
		 replace zip3`pre'=0 if zipct`pre'<1250
		 replace zip3`pre'=0 if kv_`pre'==0 
		 replace zip3`pre'=STATE_encode+1000 if zip3`pre'==0 
		 
		 capture rename mxicd8_chronic mxicd8_circulatory

		 logit `var'  i.postmerger#i.sex_pat i.vmerger##i.sex_pat i.vmerger#c.nyear i.p_age i.agecat##i.sex_pat i.qtr [pweight=samp_weight] if flag`i'==1 & kv_`pre'==1 , cluster( propScoreId)
		 
		   matrix sesq = vecdiag(e(V))' 
			svmat sesq
		  getColnames
		  getRegStats

		  rename rownames rn`pre'_`i'
		  capture drop j 
		  capture drop i
		  rename y1 beta_`pre'_`i'
		  rename sesq1 sesq_`pre'_`i'

		  nlcom exp(_cons+ _b[1.postmerger#1.sex_pat]+_b[1.vmerger] + _b[77.p_age])/(1+exp(_cons+ _b[1.postmerger#1.sex_pat]+_b[1.vmerger] + _b[77.p_age])) - (exp(_cons+  _b[77.p_age]+_b[1.vmerger])/(1+exp( _cons+ _b[1.vmerger] +_b[77.p_age])))
		   mat A=r(b)
		   mat B=r(V)
		   mat C=[A,B]
		  nlcom exp(_cons+ _b[1.postmerger#2.sex_pat]+_b[1.vmerger]+_b[1.vmerger#2.sex_pat] + _b[77.p_age] + _b[2.sex_pat]+ _b[3.agecat#2.sex_pat])/(1+exp(_cons+ _b[1.postmerger#2.sex_pat]+_b[1.vmerger]+_b[1.vmerger#2.sex_pat] + _b[2.sex_pat] + _b[77.p_age] + _b[3.agecat#2.sex_pat])) -   (exp( _cons+ _b[1.vmerger]+_b[1.vmerger#2.sex_pat]+_b[77.p_age] + _b[2.sex_pat] + _b[3.agecat#2.sex_pat])/(1+exp(_cons+ _b[1.vmerger]+_b[1.vmerger#2.sex_pat] + _b[2.sex_pat]+_b[77.p_age] + _b[3.agecat#2.sex_pat]))) 
		   mat A=r(b)
		   mat B=r(V)
		   mat D=[A,B]
		   mat E=[D',C']
		  svmat E
		  rename E1 mef_`pre'_`i'
		  rename E2 mem_`pre'_`i'

			 }
		  }
		}

save "./Overall_hyp_lim.dta", replace

	*replace state=0 if state==7 | state==39
	rename mxicd8_circulatory mxicd8_chronic 

		foreach var of varlist oc_death {
		 forval i=1(1)3  {


		if regexm("`var'","oc_([a-z0-9]+)$") { 
		  local pre = "`=regexs(1)'"

		  local j=`i'-1
		  
		 capture drop zip3`pre'`j' 
		 capture drop zip3`pre'`i' 
		 capture drop  zipct`pre'`i'
		 gen zip3`pre'`i'=zip3dig 
		 egen zipct`pre'`i'=count(zip3dig) if kv_`pre'==1 , by(zip3dig)
		 replace zip3`pre'`i'=0 if zipct`pre'`i'<1250
		 replace zip3`pre'`i'=0 if kv_`pre'==0 
		 replace zip3`pre'`i'=STATE_encode+1000 if zip3`pre'`i'==0 
		  
		 logit `var'  i.postmerger#i.sex_pat i.vmerger##i.sex_pat i.vmerger#c.nyear i.p_age i.agecat##i.sex_pat i.qtr [pweight=samp_weight] if flag`i'==1 & kv_`pre'==1 , cluster( propScoreId)
		 
		   matrix sesq = vecdiag(e(V))' 
			svmat sesq
		  getColnames
		  getRegStats

		  rename rownames rn`pre'_`i'
		  capture drop j 
		  capture drop i
		  rename y1 beta_`pre'_`i'
		  rename sesq1 sesq_`pre'_`i'

		  nlcom exp(_cons + _b[1.postmerger#1.sex_pat]+_b[1.vmerger] + _b[3.agecat])/(1+exp(_cons+_b[1.postmerger#1.sex_pat]+_b[1.vmerger] + _b[3.agecat])) - (exp(_cons+ _b[3.agecat]+ _b[1.vmerger])/(1+exp( _cons+_b[1.vmerger] +_b[3.agecat])))
		   mat A=r(b)
		   mat B=r(V)
		   mat C=[A,B]
		  nlcom exp(_cons+_b[1.postmerger#2.sex_pat]+_b[1.vmerger]+_b[1.vmerger#2.sex_pat] + _b[3.agecat] + _b[2.sex_pat]+ _b[3.agecat#2.sex_pat])/(1+exp(_cons+_b[1.postmerger#2.sex_pat]+_b[1.vmerger]+_b[1.vmerger#2.sex_pat] + _b[3.agecat] + _b[2.sex_pat] + _b[3.agecat#2.sex_pat])) -   (exp(_cons+ _b[1.vmerger]+_b[1.vmerger#2.sex_pat] + _b[3.agecat] + _b[2.sex_pat] + _b[3.agecat#2.sex_pat])/(1+exp(_cons+_b[1.vmerger]+_b[1.vmerger#2.sex_pat] + _b[3.agecat]+ _b[2.sex_pat]+ _b[3.agecat#2.sex_pat]))) 
		   mat A=r(b)
		   mat B=r(V)
		   mat D=[A,B]
		   mat E=[D',C']
		  svmat E
		  rename E1 mef_`pre'_`i'
		  rename E2 mem_`pre'_`i'
		 
		   }
		  }
		}

save "./Overall_death_lim.dta", replace
clear
 